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1. INTRODUCTION 

The representation of a dynamic system by a nonlinear stochastic model makes it possible to take 
into account factors caused by nonlinear laws of nature and perform a better analysis. To obtain a model with 
good predictive properties, informative measurement data and a suitable model structure capable of 
accurately describing the dynamics of the process are required; therefore, when constructing models of 
nonlinear systems, parametric identification methods are used. Traditionally, the maximum likelihood (ML) 
method is used to solve the problem of parametric identification [1]-[3]. In case of using dynamic models 
with Gaussian noise, the corresponding identification criterion is written on the basis of the equations of the 
extended Kalman filter (EKF) [4], [5]. Although the EKF is widely used, this filter has some drawbacks. EKF 
applies the standard technique to linearize a nonlinear model. It requires the sufficient differentiability of the 
dynamic state and the susceptibility to biasing and to divergence of the state estimates. This method is 
sub-optimal and can easily lead to the divergence. The EKF achieves only the first-order accuracy and 
produces a good result only if the initial estimation error and disturbing noises are small. 

These difficulties can be successfully overcome with such nonlinear filters as the cubature Kalman 
filter [6], [7] and the unscented Kalman filter (UKF) [8]-[12]. Julier et al. [8], [9] proposed derivative free 
alternative to the EKF for state estimation the unscented Kalman filter (UKF). The UKF has been developed 
for the case of highly nonlinear state estimation problems. The UKF performs a Gaussian approximation with 
a limited number of points (sigma points), using the unscented transform. This technique is used to linearize 
a function a random variable via the linear regression based on the points drawn from the prior distribution of 
the random variable. The UKF has the same computational complexity as the EKF has [12], but UKF does 
not require the Jacobians computing and can achieve the second-order accuracy of the Taylor expansion. 

When solving practical problems, statistical parameters of noise are set inaccurately or they are 
completely unknown. The presence of outliers in the measurement data makes the further determination of 
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such characteristics complicated. When using the incorrect a priori information about the noise properties of 
the system and/or the measurements, the obtained estimates may be biased. Covariance matrices are determined 
according to the results of the analysis of the source data or modeling. Note that the correct determination of the 
covariance matrices of noise processes affects the accuracy of the estimation of the state vector. 

A possible solution to the problem is seen in the use of adaptive methods of measuring data 
processing [13]-[17]. Adaptive algorithms allow us to jointly evaluate the state vector and covariance noise 
matrices. In this research, the sub-optimal Sage-Husa estimator [13], [16], [17] is combined with the UKF 
algorithm in order to estimate and improve the statistical properties of the process noise. Such improvement 
reduces the model error, suppresses the filtering divergence and improves the filtering accuracy. For the more 
accurate construction of a mathematical model, it is proposed to additionally evaluate the disturbances from 
the measurements of the residual differences, using a two-stage parametric identification procedure. 

This paper has the following structure. Section 2 provides a mathematical description of the motion 
model of the navigation satellite. A two-step parametric identification procedure is described in section 3. 
The results of applying the two-stage parametric identification procedure in constructing a satellite motion 
model are given in section 4. The found parameters of the solar radiation model are given. A comparison of 
the accuracy of the constructed motion models of navigation satellite based on the ML method and the 
two-stage identification procedure is also given in section 4. The conclusion is provided in section 5. 


2. MOTION MODEL OF THE NAVIGATION SATELLITE 

The quality of the ephemeris-temporal support for Global Navigation Satellite System (GNSS) 
technologies depends on adequacy of the applied mathematical models describing the orbital motion of 
navigation satellites. Consider the following stochastic nonlinear continuous-discrete model of orbital motion 
of the navigation satellite [18], [19] is (1) and (2). 


rO) = — Paar +g (rO) + g(r) + gr, 0,0) + wO,t E [tot] (1) 


f(RO) 
S(tk+1) = Tear) + Vtr), k = 0,1,...,N—1 (2) 
h(R(tpk+1)) 
r(t) 
r(t) 


T 
inertial coordinate system, r(t) = (v.(0), v(t), v,(t)) is the velocity vector of the navigation satellite in an 


Where, R(t) = ( )r@ = (x (t), y(t), z(t)’ is the coordinate vector of the navigation satellite in an 


inertial coordinate system; f(R(t)), h(R(ty41)) are functions, where u is the gravitational constant, Mg is 


the mass of the Earth; ||r(t)]|| = fx? +y) + z(t)is the radius of the orbit, g:(r(t)) is the 
perturbations, caused by the non-sphericity of the Earth's geopotential, g2 (r(t)) is the perturbations, caused 
by the gravitational influence of the Moon, the Sun and/or the other planets, g3(r(t),r(t), 0) is perturbations 
from the solar radiation (SR); 0 € Ng is the vector of unknown parameters; s(t;,41) is the measurement 
vector (for example, pseudo range, query range, satellite laser ranging (SLR) from ground points to the 
navigation spacecraft). In a particular case, a posteriori ephemeris of navigation spacecraft obtained by 
various processing centers can act as measurements (i.e. A(R (tk+1)) = r(tķ+1). In (1) and (2): 

— The random vectors w(t) and v(tęķ+1). form white Gaussian noises with unknown noises covariance 

matrices 


E[w(t)] = 0, E[w(t)w" (z)] = Qu Ct) — T) 
E[ v(te41)] = 0, E| VCtr+1) V7 tind] = Q, (tee DOK1 
E[v(tes1)w? (1)] = 0,k,i = 0,1,...,N — 1,7 E [to ty] 


— The state vector R(t)in moment todefined by 
E[R(to)] = R(to), E [(R (to) — R(to))(R(to) — RCto)) | = PC) 


and has no correlation with w(t), v(ty41). 
A description of each of the forces affecting on a satellite can be found in, for example, [19], [20]. It 
is important to note that some of these force models include parameters which numerical values are only 
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partially known. In the formation of model (1) and (2) it remains problematic to take into account 
perturbations from solar radiation pressure on the satellite [21]-[25]. To compute g3(r(t),r(t),@) in an 
inertial coordinate system, the following SR model has been used [19], [23]—[25]: 


g3(r(t), F(t), 0) = A(t) - p- 7 (r(t)) - [in - (01 + 82 cos o Cret), +C) + 63 sino (r(t), 7 (1) + 
+i,.(0, + 8; cos o (r(t), 7(t)) + 6 sin o (r(t), F(t))) 
+ i; - (0, + 8 cosa (r(t), 7(t)) + 0s sina (r(t), t) )| (3) 


Here A(r(t)) is the eclipse factor, p(r(t)) is the distance between the satellite and the Sun, a(r(t),r(t)) is 
the argument of the latitude for the navigation satellite, i, = “ee is ort in the direction of solar 
D- 


iixr(t) . ; : EON 
axr is ort normal to the Sun-satellite-Earth, iz = i, X i, 


llii xr Æl 
is ort that complements the system to the right triple of vectors. 


radiation, ||-|| is the Euclidean vector norm, i, = 


3. RESEARCH METHOD 

Informative measurement data and correctly defined model parameters affect the predictive 
properties of the model and the ability to describe the dynamics of the process. Usually, the construction of 
the model (1) consists in finding estimates of the unknown parameters of the SR model. Estimation of 
unknown parameters of model (1) and (2) is carried out according to measurement data and the selected 
identification method. The a priori assumptions allow using the ML method for the parameters estimation. 
ML estimates have practically important properties, as asymptotic efficiency, asymptotic unbiasedness, 
asymptotic normality and consistency. 

We offer a two-step procedure for parametric identification of the model (1) and (2). At the first 
stage of the procedure, the parameters of the radiation pressure model are estimated using the ML method 
based on the adaptive unscented Kalman filter. At the second stage, the parameters of the unaccounted 
perturbations model are estimated based on the results of measurements of residual differences. 

Stage 1. Building the matching model (parameter identification) 
1. Solve the problem of parametric identification based on the ML method 


=~ + aR 1oy- = 
8 =arg beds 5 kao Indet Py (tk+1) + 5 2k20 €(te41)" Py "(tha1) E(tk+1) (4) 


where €(t,41,) and Py(t,+1) are defined based on the corresponding equations of the following adaptive 
UKF [13], [16], [17]: 

Initialization: 

— Set the values 


č = 0.001,n = 2,9 =k = 0,p = 0.998 
— Define the initial values 


R (tolto) = R(to), P (tolto) = P(to), Qy (to), Q, (ty) 


— Calculate 


l = (n+) —n(n = 6 is the size of R(t)), a = + 
1 4: noes 

Bo = miia- tt T Fee B,i=1,..,2n 

a = [ao Gy, Con] 


T 
A= (: —[al... a) diag (obi wy Bon) (: —[al... al) 


2n+1 2n+1 


Fork =0,N—-—1 
Prediction: 
— Define, R(ty+1|t,), P(te+1|t,) by solving of differential (5) and (6) integration 
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qe Elto = Rs (tlt, JAR (t|tk) + Re (tlt JARs (t]tx) 
+Qw(t), QO) =Qy (teat), tk StS tear, (6) 


where the transformed set of vectors is identified as: 


Re (tlt) = [f (REEI CRIEI). If Ron (Ete) nxen) 


sigma points R#(t|t;,),i = 1,n are computed in accordance with the formula (7). 


R(t|t,), i= 0, 
R? (t|t,) = $ R(t|t,) + Vn + 1D, (t|t,), i= 1,n, (7) 


R(t|t,.) — Vn + ID,_n(t|t,), i=n+1,2n, 
Rs (t|tx) = [Ro (tlt, IR? (t|ti)I. me IRon (tlt) Inxc2n+1): 


D; is the i — th row of the lower triangular matrix obtained by the Cholesky decomposition of P(€|t,). 
Updating: 

— Find Ro(ty41|t,) using (7) with the substitution t = ty41. 

— Calculate 


Sn(teralte) = [A(R tk+1 lt) [ACRE (tested) [--- [ACRen (tesa ti) I nsccomesy’ 
ettigi) = S(tk+1) — Sn (tr+1ltx)a, 
-p 


ESE 
Ôn Cte) = (1 — %) (ty) + TelEltk+1)ET(tk+1) — z 
= Xito Éi (h (RÈ (tesalte)) = Su(terrlteda) (h (RE (tesalti)) = Su(terrlte)a) } 


Ps (tyes) = Sp(tesrlti DASA Ctk+1ltk) + Qy (tera, 
Prs(tk+1) = RsCtr+1ltr) AS (tkyilti). 


— State Å(tk+1ltk+1) and covariance estimates P(t,+1|t,41) are computed according to the equations 


K(tk+1) a Pps (thsi )P5* (ters) 
R(tysiltess) = ROesilte) + K (tk+1)Eltk+1) 
P(tk+1ltk+1) = P(tk+1ltr) = K (tk41)Ps(tk+41)K" tear) 


Ow (nea) = (= Tr) Qw (th) + fex [ctete CRT Ca + P(tk+1ltk+1) — 


2 pi (F (RE Ctusrlte)) — R ltrs lte) ) (F (RF Celt) — Ress i) | 
End. 


Note that the cost function in (4) is known to have many local optima. Often to solve the 
optimization problem (4) are used Newton’s method and various quasi-Newton methods, which are the local 
ones. These methods are sensitive to the setting of initial values and the accuracy of determining the gradient. 
In general, it is impossible to determine what influenced the low accuracy of the model, the inaccurate 
definition of the model structure or the found local minimum. A reasonable approach to optimization 
problem is to use global optimization methods. In the paper for the solution (4) is used global optimization 
approach based on the sequential quadratic programming method. 

Stage 2. Specifying the matching model (identification of unaccounted disturbances from measurements of 
residual differences) 
— Calculate the residual differences A(t,1) based on the estimates 6 obtained in step 1. 


A( titi) = S(tk4+1) = $1 (teas) 


Two-stage parametric identification procedure to predict satellite ... (Oksana Sergeevna Chernikova) 


5352. 0O ISSN: 2088-8708 


where 
§4(tear) = h(R(tesaltess)) (8) 


R(ty+1\th41)- the estimate of the state vector obtained at the first stage. 
— Choose the following model 


Ai (thsi) = ab + ality, + abt2,, + bÍ cos (m=) + 
bésin (=) + cicos (=) + cÍ sin (=) (9) 
i = 1,...,m, m - the size of s(t,41) (here m=3), T-defined value, and estimate of the unknown parameters 


ő = (ab, al, aż, bi, bi, gi ch) ,i = 1,2,3, using least squares method [26]: 


6 = (BB) *BTA(ty41), 


2 2mty 8 2mty Ant, a Ant, 
1 ty tí cos T sın T cos m sın T 


where B = [e > os = a va 
1 ty teens (=) Si (=x) cos (sin (=) 


— Calculate Â(tp41) taking into account the estimates 6 found by the (9). 
— Calculate 


$5 (thea) = Si (Gia) + AG) k = 0, ..,N — 1 (10) 


using A(t,41) and ŝ; (tk+1) found by the (8). 


4. RESULTS AND DISCUSSION 

As the measurement data we have taken the rapid ephemeris of the GPS from 07/14/2016, obtained 
by the international GNSS service. In this case, the satellite makes more than one revolution around the Earth 
(passes through the different light zones). At the initial time, we compute the velocity of the satellite on the 
basis of rapid ephemeris using Everett interpolation. Estimation of the SR parameters (3) can be carried out 
using the ML method according to the trajectory observations in areas of total illumination and penumbra 
zones [27]. The quality of the parameters estimates found was determined by the accuracy of predicting the 
navigation satellite motion trajectory: 


RMSE} = [EENES Cte) — SiC) 
$ 1 S ¥ ne 2 
RMSE; = [A zys (thai) — $? Cte) 


Here {s(ty41),k = 0,1,...,N — 1} is final ephemeris received on July 14, 2016, {8} (tk+1), k = 


0,1,...,N — 1} is the predicted trajectory on 06/14/2016 for the filter equation at @ and {$3(t,41),k = 
0,1,..., N — 1} is the obtained trajectory on 06/14/2016 based on the found parameters é, S (tga) is final 
ephemeris received on 06/15/2016, {ŝi (tk+1), k = 0,1,..., N — 1} is the obtained trajectory on 06/15/2016 
for the filter equation at Ê and {$3(t,41),k = 0,1,..., N — 1}is the obtained trajectory on 06/15/2016 based 


on the found parameters 9. The obtained results are presented in Table 1. 

Thus, the result of the sunlight SR parameters specification is that it is possible to increase the 
accuracy of the satellite trajectory prediction. The results obtained show that the model based on the two- 
stage parametric identification procedure is more accurate than the model obtained using the one-stage 
parametric identification procedure. 
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Table 1. Results of estimating parameters of solar radiation model 


One-stage parametric Two-stage parametric identification procedure 
identification procedure 
1.06906362 0.0955 0.0668 —0.0143 
0.06858372 —0.0010 —0.0051 0.0007 
0.04729392 7 0.0000 0.0001 0.0000 
_ _|0.11131100 ð = 1077 x]—0.2458 —0.1025 0.1188 
Parameters estimates ð = | 0.06708731 0.1703 —0.1729 0.1879 
0.09272575 —0.0009 —0.0051 0.0061 
0.09483054 —0.0186 0.0114 —0.0141 
0.13523427 
0.10924206 
Evaluation of the quality of building a RMSE} = 3.2701 x 10-8 RMSE? = 8.1142 x 10-° 
mathematical model (km) 
Evaluation of the prediction quality (km) RMSE} = 3.3985 x 1078 RMSE = 1.9186 x 1078 


5. CONCLUSION 

In this paper the statistic estimator based on the adaptive modification UKF with noise is used for 
the SR model parameters estimation. This approach improves the robustness of the traditional UKF with the 
respect to the variable noise distribution. The results show that in case of time-varying or uncertain noise 
characteristics the adaptive UKF is more efficient than the conventional UKF in terms of the fast 
convergence and the state estimation accuracy. 

Offered approach to the construction of a satellite motion model and the prediction of ephemeris is 
based on a complex application of the maximum likelihood method, a nonlinear filtering algorithm and the 
identification of a complex component of the satellite motion model. The use of a two-stage parametric 
identification procedure that combines the estimation of the parameters of the matching SR model from 
ephemerides with the refinement of the residual accelerations made it possible to increase the accuracy of 
determining the unknown parameters of the satellite motion model. In the future, to improve the accuracy of 
the model construction, a two-stage gradient identification procedure based on a combination of several 
adaptive filters will be developed. 
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